#ifndef __CTP33D3D__
#define __CTP33D3D__
#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// P3 element, conforming, 3D
// ***********************************************************************
static INT C_T_P3_3D_3D_dof[4] = {3,6,3,0 };
static INT C_T_P3_3D_3D_Num_Bas = 60;
static INT C_T_P3_3D_3D_Value_Dim = 3;
static INT C_T_P3_3D_3D_Inter_Dim = 3;
static INT C_T_P3_3D_3D_Polydeg = 3;  
static bool C_T_P3_3D_3D_IsOnlyDependOnRefCoord = 1;
static INT C_T_P3_3D_3D_Accuracy = 3; 
static MAPTYPE C_T_P3_3D_3D_Maptype = Affine;

static void C_T_P3_3D_3D_InterFun(DOUBLE RefCoord[ ], INT dim, DOUBLE *values)
{
    DOUBLE x, y, z;
    x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];

    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0*x*y - (11.0*y)/2.0 - (11.0*z)/2.0 - (11.0*x)/2.0 + 18.0*x*z + 18.0*y*z - (27.0*x*y*y)/2.0 
                    - (27.0*x*x*y)/2.0 - (27.0*x*z*z)/2.0 - (27.0*x*x*z)/2.0 - (27.0*y*z*z)/2.0 - (27.0*y*y*z)/2.0 
                    + 9.0*x*x - (9.0*x*x*x)/2.0 + 9.0*y*y - (9.0*y*y*y)/2.0 + 9.0*z*z - (9.0*z*z*z)/2.0 - 27.0*x*y*z + 1.0;
        values[9*1+i*4] = x - (9.0*x*x)/2.0 + (9.0*x*x*x)/2.0;
        values[9*2+i*4] = y - (9.0*y*y)/2.0 + (9.0*y*y*y)/2.0;
        values[9*3+i*4] = z - (9.0*z*z)/2.0 + (9.0*z*z*z)/2.0;
        values[9*4+i*4] = 9.0*x - (45.0*x*y)/2.0 - (45.0*x*z)/2.0 + (27.0*x*y*y)/2.0 + 27.0*x*x*y + (27.0*x*z*z)/2.0 + 27.0*x*x*z 
                    - (45.0*x*x)/2.0 + (27.0*x*x*x)/2.0 + 27.0*x*y*z;
        values[9*5+i*4] = (9.0*x*y)/2.0 - (9.0*x)/2.0 + (9.0*x*z)/2.0 - (27.0*x*x*y)/2.0 - (27.0*x*x*z)/2.0 + 18.0*x*x - (27.0*x*x*x)/2.0;
        values[9*6+i*4] = 9.0*y - (45.0*x*y)/2.0 - (45.0*y*z)/2.0 + 27.0*x*y*y + (27.0*x*x*y)/2.0 + (27.0*y*z*z)/2.0 + 27.0*y*y*z 
                    - (45.0*y*y)/2.0 + (27.0*y*y*y)/2.0 + 27.0*x*y*z;
        values[9*7+i*4] = (9.0*x*y)/2.0 - (9.0*y)/2.0 + (9.0*y*z)/2.0 - (27.0*x*y*y)/2.0 - (27.0*y*y*z)/2.0 + 18.0*y*y - (27.0*y*y*y)/2.0;
        values[9*8+i*4] = 9.0*z - (45.0*x*z)/2.0 - (45.0*y*z)/2.0 + 27.0*x*z*z + (27.0*x*x*z)/2.0 + 27.0*y*z*z + (27.0*y*y*z)/2.0 
                    - (45.0*z*z)/2.0 + (27.0*z*z*z)/2.0 + 27.0*x*y*z;
        values[9*9+i*4] = (9.0*x*z)/2.0 - (9.0*z)/2.0 + (9.0*y*z)/2.0 - (27.0*x*z*z)/2.0 - (27.0*y*z*z)/2.0 + 18.0*z*z - (27.0*z*z*z)/2.0;
        values[9*10+i*4] = (27.0*x*x*y)/2.0 - (9.0*x*y)/2.0;
        values[9*11+i*4] = (27.0*x*y*y)/2.0 - (9.0*x*y)/2.0;
        values[9*12+i*4] = (27.0*x*x*z)/2.0 - (9.0*x*z)/2.0;
        values[9*13+i*4] = (27.0*x*z*z)/2.0 - (9.0*x*z)/2.0;
        values[9*14+i*4] = (27.0*y*y*z)/2.0 - (9.0*y*z)/2.0;    
        values[9*15+i*4] = (27.0*y*z*z)/2.0 - (9.0*y*z)/2.0;
        values[9*16+i*4] = 27.0*x*y*z;
        values[9*17+i*4] = 27.0*y*z - 27.0*y*z*z - 27.0*y*y*z - 27.0*x*y*z;
        values[9*18+i*4] = 27.0*x*z - 27.0*x*z*z - 27.0*x*x*z - 27.0*x*y*z;
        values[9*19+i*4] = 27.0*x*y - 27.0*x*y*y - 27.0*x*x*y - 27.0*x*y*z;
    }
}

static void C_T_P3_3D_3D_D000(ELEMENT *Elem, DOUBLE Coord[ ], DOUBLE RefCoord[ ], DOUBLE *values)
{
    DOUBLE x, y, z;
    x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    //x =1.0/3.0; y=1.0/3.0; z=0.0;
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0*x*y - (11.0*y)/2.0 - (11.0*z)/2.0 - (11.0*x)/2.0 + 18.0*x*z + 18.0*y*z - (27.0*x*y*y)/2.0 
                    - (27.0*x*x*y)/2.0 - (27.0*x*z*z)/2.0 - (27.0*x*x*z)/2.0 - (27.0*y*z*z)/2.0 - (27.0*y*y*z)/2.0 
                    + 9.0*x*x - (9.0*x*x*x)/2.0 + 9.0*y*y - (9.0*y*y*y)/2.0 + 9.0*z*z - (9.0*z*z*z)/2.0 - 27.0*x*y*z + 1.0;
        values[9*1+i*4] = x - (9.0*x*x)/2.0 + (9.0*x*x*x)/2.0;
        values[9*2+i*4] = y - (9.0*y*y)/2.0 + (9.0*y*y*y)/2.0;
        values[9*3+i*4] = z - (9.0*z*z)/2.0 + (9.0*z*z*z)/2.0;
        values[9*4+i*4] = 9.0*x - (45.0*x*y)/2.0 - (45.0*x*z)/2.0 + (27.0*x*y*y)/2.0 + 27.0*x*x*y + (27.0*x*z*z)/2.0 + 27.0*x*x*z 
                    - (45.0*x*x)/2.0 + (27.0*x*x*x)/2.0 + 27.0*x*y*z;
        values[9*5+i*4] = (9.0*x*y)/2.0 - (9.0*x)/2.0 + (9.0*x*z)/2.0 - (27.0*x*x*y)/2.0 - (27.0*x*x*z)/2.0 + 18.0*x*x - (27.0*x*x*x)/2.0;
        values[9*6+i*4] = 9.0*y - (45.0*x*y)/2.0 - (45.0*y*z)/2.0 + 27.0*x*y*y + (27.0*x*x*y)/2.0 + (27.0*y*z*z)/2.0 + 27.0*y*y*z 
                    - (45.0*y*y)/2.0 + (27.0*y*y*y)/2.0 + 27.0*x*y*z;
        values[9*7+i*4] = (9.0*x*y)/2.0 - (9.0*y)/2.0 + (9.0*y*z)/2.0 - (27.0*x*y*y)/2.0 - (27.0*y*y*z)/2.0 + 18.0*y*y - (27.0*y*y*y)/2.0;
        values[9*8+i*4] = 9.0*z - (45.0*x*z)/2.0 - (45.0*y*z)/2.0 + 27.0*x*z*z + (27.0*x*x*z)/2.0 + 27.0*y*z*z + (27.0*y*y*z)/2.0 
                    - (45.0*z*z)/2.0 + (27.0*z*z*z)/2.0 + 27.0*x*y*z;
        values[9*9+i*4] = (9.0*x*z)/2.0 - (9.0*z)/2.0 + (9.0*y*z)/2.0 - (27.0*x*z*z)/2.0 - (27.0*y*z*z)/2.0 + 18.0*z*z - (27.0*z*z*z)/2.0;
        values[9*10+i*4] = (27.0*x*x*y)/2.0 - (9.0*x*y)/2.0;
        values[9*11+i*4] = (27.0*x*y*y)/2.0 - (9.0*x*y)/2.0;
        values[9*12+i*4] = (27.0*x*x*z)/2.0 - (9.0*x*z)/2.0;
        values[9*13+i*4] = (27.0*x*z*z)/2.0 - (9.0*x*z)/2.0;
        values[9*14+i*4] = (27.0*y*y*z)/2.0 - (9.0*y*z)/2.0;    
        values[9*15+i*4] = (27.0*y*z*z)/2.0 - (9.0*y*z)/2.0;
        values[9*16+i*4] = 27.0*x*y*z;
        values[9*17+i*4] = 27.0*y*z - 27.0*y*z*z - 27.0*y*y*z - 27.0*x*y*z;
        values[9*18+i*4] = 27.0*x*z - 27.0*x*z*z - 27.0*x*x*z - 27.0*x*y*z;
        values[9*19+i*4] = 27.0*x*y - 27.0*x*y*y - 27.0*x*x*y - 27.0*x*y*z;
    }
}

static void C_T_P3_3D_3D_D100(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0*x + 18.0*y + 18.0*z - 27.0*x*y - 27.0*x*z - 27.0*y*z - (27.0*x*x)/2.0 - (27.0*y*y)/2.0 - (27.0*z*z)/2.0 - 11.0/2.0;
        values[9*1+i*4] = (27.0*x*x)/2.0 - 9.0*x + 1.0;
        values[9*2+i*4] = 0.0;
        values[9*3+i*4] = 0.0;
        values[9*4+i*4] = 54.0*x*y - (45.0*y)/2.0 - (45.0*z)/2.0 - 45.0*x + 54.0*x*z + 27.0*y*z + (81.0*x*x)/2.0 + (27.0*y*y)/2.0 + (27.0*z*z)/2.0 + 9.0;
        values[9*5+i*4] = 36.0*x + (9.0*y)/2.0 + (9.0*z)/2.0 - 27.0*x*y - 27.0*x*z - (81.0*x*x)/2.0 - 9.0/2.0;
        values[9*6+i*4] = 27.0*x*y - (45.0*y)/2.0 + 27.0*y*z + 27.0*y*y;
        values[9*7+i*4] = (9.0*y)/2.0 - (27.0*y*y)/2.0;
        values[9*8+i*4] = 27.0*x*z - (45.0*z)/2.0 + 27.0*y*z + 27.0*z*z;
        values[9*9+i*4] = (9.0*z)/2.0 - (27.0*z*z)/2.0;
        values[9*10+i*4] = 27.0*x*y - (9.0*y)/2.0;
        values[9*11+i*4] = (27.0*y*y)/2.0 - (9.0*y)/2.0;
        values[9*12+i*4] = 27.0*x*z - (9.0*z)/2.0;
        values[9*13+i*4] = (27.0*z*z)/2.0 - (9.0*z)/2.0;
        values[9*14+i*4] = 0.0;
        values[9*15+i*4] = 0.0;
        values[9*16+i*4] = 27.0*y*z;
        values[9*17+i*4] = -27.0*y*z;
        values[9*18+i*4] = 27.0*z - 54.0*x*z - 27.0*y*z - 27.0*z*z;
        values[9*19+i*4] = 27.0*y - 54.0*x*y - 27.0*y*z - 27.0*y*y;
    }
}

static void C_T_P3_3D_3D_D010(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0*x + 18.0*y + 18.0*z - 27*x*y - 27*x*z - 27*y*z - (27*x*x)/2.0 - (27*y*y)/2.0 - (27*z*z)/2.0 - 11.0/2.0;
        values[9*1+i*4] = 0.0;
        values[9*2+i*4] = (27.0*y*y)/2.0 - 9.0*y + 1.0;
        values[9*3+i*4] = 0.0;
        values[9*4+i*4] = 27.0*x*y - (45*x)/2.0 + 27*x*z + 27*x*x;
        values[9*5+i*4] = (9.0*x)/2.0 - (27.0*x*x)/2.0;
        values[9*6+i*4] = 54.0*x*y - 45*y - (45*z)/2.0 - (45*x)/2.0 + 27*x*z + 54*y*z + (27*x*x)/2.0 + (81*y*y)/2.0 + (27*z*z)/2.0 + 9;
        values[9*7+i*4] = (9.0*x)/2.0 + 36.0*y + (9.0*z)/2.0 - 27.0*x*y - 27.0*y*z - (81.0*y*y)/2.0 - 9.0/2.0;
        values[9*8+i*4] = 27.0*x*z - (45*z)/2.0 + 27*y*z + 27*z*z;
        values[9*9+i*4] = (9.0*z)/2.0 - (27.0*z*z)/2.0;
        values[9*10+i*4] = (27.0*x*x)/2.0 - (9*x)/2.0;
        values[9*11+i*4] = 27.0*x*y - (9.0*x)/2.0;
        values[9*12+i*4] = 0.0;
        values[9*13+i*4] = 0.0;
        values[9*14+i*4] = 27.0*y*z - (9.0*z)/2.0;    
        values[9*15+i*4] = (27.0*z*z)/2.0 - (9.0*z)/2.0;
        values[9*16+i*4] = 27.0*x*z;
        values[9*17+i*4] = 27.0*z - 27.0*x*z - 54.0*y*z - 27.0*z*z;
        values[9*18+i*4] = -27.0*x*z;
        values[9*19+i*4] = 27.0*x - 54.0*x*y - 27.0*x*z - 27.0*x*x;
    }
}
static void C_T_P3_3D_3D_D001(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.*x + 18.*y + 18.*z - 27.*x*y - 27.*x*z - 27.*y*z - (27.*x*x)/2.0 - (27.*y*y)/2.0 - (27.*z*z)/2.0 - 11.0/2.0;
        values[9*1+i*4] = 0;
        values[9*2+i*4] = 0;
        values[9*3+i*4] = (27.0*z*z)/2.0 - 9*z + 1.0;
        values[9*4+i*4] = 27.0*x*y - (45.0*x)/2.0 + 27.0*x*z + 27.0*x*x;
        values[9*5+i*4] = (9*x)/2.0 - (27*x*x)/2.0;
        values[9*6+i*4] = 27.0*x*y - (45.0*y)/2.0 + 27.0*y*z + 27.0*y*y;
        values[9*7+i*4] = (9*y)/2.0 - (27*y*y)/2.0;
        values[9*8+i*4] = 27.0*x*y - (45.0*y)/2.0 - 45.0*z - (45.0*x)/2.0 + 54.0*x*z + 54.0*y*z + (27.0*x*x)/2.0 + (27.0*y*y)/2.0 + (81.0*z*z)/2.0 + 9.0;
        values[9*9+i*4] = (9.0*x)/2.0 + (9.0*y)/2.0 + 36.0*z - 27.0*x*z - 27.0*y*z - (81.0*z*z)/2.0 - 9.0/2.0;
        values[9*10+i*4] = 0.0;
        values[9*11+i*4] = 0.0;
        values[9*12+i*4] = (27.0*x*x)/2.0 - (9.0*x)/2.0;
        values[9*13+i*4] = 27.0*x*z - (9.0*x)/2.0;
        values[9*14+i*4] = (27.0*y*y)/2.0 - (9*y)/2.0;   
        values[9*15+i*4] = 27.0*y*z - (9.0*y)/2.0;
        values[9*16+i*4] = 27.0*x*y;
        values[9*17+i*4] = 27.0*y - 27.0*x*y - 54.0*y*z - 27.0*y*y;
        values[9*18+i*4] = 27.0*x - 27.0*x*y - 54.0*x*z - 27.0*x*x;
        values[9*19+i*4] = -27.0*x*y;
    }
}
static void C_T_P3_3D_3D_D200(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0 - 27.0*y - 27.0*z - 27.0*x;
        values[9*1+i*4] = 27.0*x - 9.0;
        values[9*2+i*4] = 0.0;
        values[9*3+i*4] = 0.0;
        values[9*4+i*4] = 81.0*x + 54.0*y + 54.0*z - 45.0;
        values[9*5+i*4] = 36.0 - 27.0*y - 27.0*z - 81.0*x;
        values[9*6+i*4] = 27.0*y;
        values[9*7+i*4] = 0.0;
        values[9*8+i*4] = 27.0*z;
        values[9*9+i*4] = 0.0;
        values[9*10+i*4] = 27.0*y;
        values[9*11+i*4] = 0.0;
        values[9*12+i*4] = 27.0*z;
        values[9*13+i*4] = 0.0;
        values[9*14+i*4] = 0.0;   
        values[9*15+i*4] = 0.0;
        values[9*16+i*4] = 0.0;
        values[9*17+i*4] = 0.0;
        values[9*18+i*4] = -54.0*z;
        values[9*19+i*4] = -54.0*y;
    }
}
static void C_T_P3_3D_3D_D020(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0 - 27.0*y - 27.0*z - 27.0*x;
        values[9*1+i*4] = 0.0;
        values[9*2+i*4] = 27.0*y - 9.0;
        values[9*3+i*4] = 0.0;
        values[9*4+i*4] = 27.0*x;
        values[9*5+i*4] = 0.0;   
        values[9*6+i*4] = 54.0*x + 81.0*y + 54*z - 45.0;
        values[9*7+i*4] = 36.0 - 81.0*y - 27.0*z - 27.0*x;
        values[9*8+i*4] = 27.0*z;
        values[9*9+i*4] = 0.0;
        values[9*10+i*4] = 0.0;
        values[9*11+i*4] = 27.0*x;
        values[9*12+i*4] = 0.0;
        values[9*13+i*4] = 0.0;
        values[9*14+i*4] = 27.0*z;   
        values[9*15+i*4] = 0.0;
        values[9*16+i*4] = 0.0;
        values[9*17+i*4] = -54.0*z;
        values[9*18+i*4] = 0.0;
        values[9*19+i*4] = -54.0*x;
    }
}
static void C_T_P3_3D_3D_D002(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0 - 27.0*y - 27.0*z - 27.0*x;
        values[9*1+i*4] = 0.0;
        values[9*2+i*4] = 0.0;
        values[9*3+i*4] = 27.0*z - 9.0;
        values[9*4+i*4] = 27.0*x;
        values[9*5+i*4] = 0.0;
        values[9*6+i*4] = 27.0*y;
        values[9*7+i*4] = 0.0;
        values[9*8+i*4] = 54.0*x + 54.0*y + 81.0*z - 45.0;
        values[9*9+i*4] = 36.0 - 27.0*y - 81.0*z - 27.0*x;
        values[9*10+i*4] = 0.0;
        values[9*11+i*4] = 0.0;
        values[9*12+i*4] = 0.0;
        values[9*13+i*4] = 27.0*x;
        values[9*14+i*4] = 0.0;  
        values[9*15+i*4] = 27.0*y;
        values[9*16+i*4] = 0.0;
        values[9*17+i*4] = -54.0*y;
        values[9*18+i*4] = -54.0*x;
        values[9*19+i*4] = 0.0;
    }
}
static void C_T_P3_3D_3D_D110(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0 - 27.0*y - 27.0*z - 27.0*x;
        values[9*1+i*4] = 0.0;
        values[9*2+i*4] = 0.0;
        values[9*3+i*4] = 0.0;
        values[9*4+i*4] = 54.0*x + 27.0*y + 27.0*z - 22.5;
        values[9*5+i*4] = 9.0/2.0 - 27.0*x;
        values[9*6+i*4] = 27.0*x + 54.0*y + 27.0*z - 22.5;
        values[9*7+i*4] = 9/2.0 - 27.0*y;
        values[9*8+i*4] = 27.0*z;
        values[9*9+i*4] = 0.0;
        values[9*10+i*4] = 27.0*x - 4.5;
        values[9*11+i*4] = 27.0*y - 4.5;
        values[9*12+i*4] = 0.0;
        values[9*13+i*4] = 0.0;
        values[9*14+i*4] = 0.0;   
        values[9*15+i*4] = 0.0;
        values[9*16+i*4] = 27.0*z;
        values[9*17+i*4] = -27.0*z;
        values[9*18+i*4] = -27.0*z;
        values[9*19+i*4] =  27.0 - 54.0*y - 27.0*z - 54.0*x;
    }
}
static void C_T_P3_3D_3D_D101(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0 - 27.0*y - 27.0*z - 27.0*x;
        values[9*1+i*4] = 0.0;
        values[9*2+i*4] = 0.0;
        values[9*3+i*4] = 0.0;
        values[9*4+i*4] = 54.0*x + 27.0*y + 27.0*z - 22.5;
        values[9*5+i*4] = 4.5 - 27.0*x;
        values[9*6+i*4] = 27.0*y;
        values[9*7+i*4] = 0.0;
        values[9*8+i*4] = 27.0*x + 27.0*y + 54.0*z - 22.5;
        values[9*9+i*4] = 4.5 - 27.0*z;
        values[9*10+i*4] = 0.0;
        values[9*11+i*4] = 0.0;
        values[9*12+i*4] = 27.0*x - 4.5;
        values[9*13+i*4] = 27.0*z - 4.5;
        values[9*14+i*4] = 0.0;   
        values[9*15+i*4] = 0.0;
        values[9*16+i*4] = 27.0*y;
        values[9*17+i*4] = -27.0*y;
        values[9*18+i*4] = 27.0 - 27.0*y - 54.0*z - 54.0*x;
        values[9*19+i*4] = -27.0*y;
    }
}
static void C_T_P3_3D_3D_D011(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    DOUBLE x, y, z;
    x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
    INT i; 
    for(i=0;i<180;i++)  values[i] = 0;
    for(i=0;i<3;i++)
    {
        values[9*0+i*4] = 18.0 - 27.0*y - 27.0*z - 27.0*x;
        values[9*1+i*4] = 0.0;
        values[9*2+i*4] = 0.0;
        values[9*3+i*4] = 0.0;
        values[9*4+i*4] = 27.0*x;
        values[9*5+i*4] = 0.0;
        values[9*6+i*4] = 27.0*x + 54.0*y + 27.0*z - 22.5;
        values[9*7+i*4] = 4.5 - 27.0*y;
        values[9*8+i*4] = 27.0*x + 27.0*y + 54.0*z - 22.5;
        values[9*9+i*4] = 4.5 - 27.0*z;
        values[9*10+i*4] = 0.0;
        values[9*11+i*4] = 0.0;
        values[9*12+i*4] = 0.0;
        values[9*13+i*4] = 0.0;
        values[9*14+i*4] = 27.0*y - 4.5;   
        values[9*15+i*4] = 27.0*z - 4.5;
        values[9*16+i*4] = 27.0*x;
        values[9*17+i*4] = 27.0 - 54.0*y - 54.0*z - 27.0*x;
        values[9*18+i*4] = -27.0*x;
        values[9*19+i*4] = -27.0*x;
    }
} 
// values of the derivatives in RefCoord[1]-RefCoord[1] direction
// values of the derivatives in RefCoord[1]-RefCoord[1] direction
static void C_T_P3_3D_3D_Nodal(ELEMENT * elem, FUNCTIONVEC *fun, INT dim, DOUBLE* values)
{
    DOUBLE coord[3];
    INT i, j;
    //先计算4个节点上的函数值
    for(i=0;i<4;i++)
    {
       // dof on the verts 
        coord[0] = elem->Vert_X[i]; 
        coord[1] = elem->Vert_Y[i];
        coord[2] = elem->Vert_Z[i];
        fun(coord, dim, values+i*dim);
    }
    INT vert0[6] = {0, 0, 0, 1, 1, 2};
    INT vert1[6] = {1, 2, 3, 2, 3, 3};
    DOUBLE W[2]  = {1.0/3.0, 2.0/3.0}; 
    //六条边上的值
     for(i=0;i<6;i++)     
     {
         for(j=0;j<2;j++)
         {
            coord[0] = (1.0-W[j])*elem->Vert_X[vert0[i]]+W[j]*elem->Vert_X[vert1[i]]; 
            coord[1] = (1.0-W[j])*elem->Vert_Y[vert0[i]]+W[j]*elem->Vert_Y[vert1[i]];
            coord[2] = (1.0-W[j])*elem->Vert_Z[vert0[i]]+W[j]*elem->Vert_Z[vert1[i]]; 
            fun(coord, dim, values+(4+i*2+j)*dim);
         }
     } 
     //四个面上的值
     DOUBLE facew[3] = {1.0/3.0, 1.0/3.0, 1.0/3.0};
     INT face0[4] = {1, 0, 0, 0};
     INT face1[4] = {2, 3, 1, 2};
     INT face2[4] = {3, 2, 3, 1};
     for(i=0;i<4;i++)
     {
        coord[0] = facew[0]*elem->Vert_X[face0[i]]+facew[1]*elem->Vert_X[face1[i]] + facew[2]*elem->Vert_X[face2[i]]; 
        coord[1] = facew[0]*elem->Vert_Y[face0[i]]+facew[1]*elem->Vert_Y[face1[i]] + facew[2]*elem->Vert_Y[face2[i]]; 
        coord[2] = facew[0]*elem->Vert_Z[face0[i]]+facew[1]*elem->Vert_Z[face1[i]] + facew[2]*elem->Vert_Z[face2[i]]; 
        fun(coord, dim, values+(16+i)*dim);
     }      
}
#endif